import numpy as np
import pandas as pd
from statsmodels.formula.api import ols
from scipy.optimize import minimize

# 假设我们有一个数据集 data.csv，包含输入变量 x1, x2 和响应变量 y
data = pd.read_csv('data1.csv',encoding='gbk')

# 构建二阶响应面模型
formula = '断裂强力N ~ x1 + x2 + x3 +I(x3**2) + I(x3**3)'
model = ols(formula, data).fit()

# 打印模型摘要
print(model.summary())

# 定义目标函数，使用二阶响应面模型的系数
def response_surface(x):
    x1, x2 ,x3 = x
    return -(model.params['Intercept'] + model.params['x1']*x1 + model.params['x2']*x2 +model.params['x3']*x3+
            model.params['I(x3 ** 2)']*x3**2 + model.params['I(x3 ** 3)']*x3**3)
# 找到响应面模型上的最优输入变量组合
initial_guess = [20, 100, 10]
bounds = [(15, 30), (100, 130), (0, 30)]
result = minimize(response_surface, initial_guess, bounds=bounds)
optimal_x = result.x
optimal_y = response_surface(optimal_x)

print(f"Optimal input variables: x1={optimal_x[0]:.4f}, x2={optimal_x[1]:.4f},x3={optimal_x[2]:.4f}")
print(f"Optimal response value: y={optimal_y:.4f}")